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ABSTRACT 

In this work, we present EnzoMHD, the extension of the cosmological code Enzo to 
include the effects magnetic fields through the ideal MHD approximation. We use a 
higher order Godunov Riemann solver for the computation of interface fluxes. We use 
two constrained transport methods to compute the electric field from those interface 
fluxes, which simultaneously advances the induction equation and maintains the diver- 
gence of the magnetic field. A third order divergence free reconstruction technique is 
used to interpolate the magnetic fields in the block structured AMR framework already 
extant in Enzo. This reconstruction also preserves the divergence of the magnetic field 
to machine precision. We use operator splitting to include gravity and cosmological ex- 
pansion. We then present a series of cosmological and non cosmological tests problems 
to demonstrate the quality of solution resulting from this combination of solvers. 



1. Introduction 



Enzo is an adaptive mesh refinement (AMR), grid-based hybrid code (hydro + N-Body) which 
is designed t o do simulations o f cosm ological structure formation. It uses the block-structured AMR 
algorithm of lBerger fc Colellal (|1989n to improve spatial resolution where required, such as in grav- 
itationally collapsing objects. The method is attractive for cosmological applications because it: 1) 
is spatially- and time-adaptive, 2) uses accurate and well-tested grid-based methods for solving the 
hydrodynamics equations and 3) can be well optimized and parallelized. The central idea behind 
AMR is to solve the evolution equations on a fixed resolution grid, adding finer grids in regions 
that require enhanced resolution. Mesh refinement can be continued to an arbitrary level, based 
on criteria involving any combination of (dark-matter and/or baryon) over density, Jeans length, 
cooling time, etc, enabling users to tailor the adaptivity to the problem of interest. Enzo solves 
the following physi cs models: collisionless dark -matter and star particles, using the particle-mesh 
N-body technique (IHockney &; Eastwood! Il985l ): gravity, using FFTs on the root grid and multi- 
grid relaxation o n the subgrids; cosm i c exp ansion; gas dynamics, using the piecewise p a raboli c 
method (PPM) (jColella fc Woodward! 11984 ) as extended to cosmology by iBryan et al.l (|1995l ): 
multi-species non-equilibrium ionization and H2 chemistry, using backward Euler time differencing 
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Anninos et al .11 19971) ; radi ative heating and cooling, using subcycled forward Euler t ime differencing 



Anninos h Norman! Il994l ): and a parameterized star formation/feedback recipe (jCen h Ostrikei 



19931 ). Enzo has been s u ccessfully used in many c osmological applications, including star for- 
mat ion (lAbel et al.1 hood I2OO2I: b'Shea et al.1 hoosl : lo'Shea fc Normanl 12007m. Lvman -aloha for- 



est (jBryan et al 



1999 



Jena et al 



galaxy clusters (jBryan fc Norman 



20051) , interstellar medium (IKritsuk h Norman 



1998 



Loken et al 



2002 



Motl et al 



2004 



2002 



Hallman et al. 



2004) and 



2006) 



More informations about Enzo are available at http://lca.ucsd.edu/projects/enzo 



One important piece of physics that is missing from this list is a proper treatment of magnetic 
fields. Magnetic fields have a broad range of impacts in a broad range of physical situations, from 
galaxy clusters to protostellar core formation. Magnetic forces can shape morphology of objects by 
forcing flow along the field lines. They can alter the energy balance by providing sources of pressure 
and energy. They can alter cooling rates by trapping electrons. Alfven waves can redistribute 
angular momentum throughout an object. They create X-ray cavities seen in some galaxy clusters. 
They accelerate cosmic rays, which play a crucial role in the energy balance of the galaxy and galaxy 
clusters. They also play a role in galactic star formation, potentially removing angular momentum 
from collapsing objects and launching protostellar winds. Creating a functional cosmological MHD 
code takes more than a single algorithm. The purpose of this paper is to document the construction 
and performance of the algorithms that will be used in MHD simula tions with Enzo in the future, 
as well as simulations that have already been done iXu et al.ll2008al ibh 



EnzoMHD is also a purpose code. In this paper, we will discuss it as a cosmological code, but 
all the same machinery applies in non-cosmological mode. All algorithms used here reduce to the 
non-cosmological limit by setting a — ► l,d — ► 0, and a — > 0. This removes any frame dependent 
terms in the equations. 

We will describe the numerical procedures in section [2j present test problems in section [3J and 
present conclusions and future plans in section [H In appendix O we present a simplified schematic 
to unify the pieces of the solver, and in appendix [A] and [B] we expand on some of the more complex 
numerical procedures. 
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2. Numerics 
2.1. Cosmological MHD Equations 

EnzoMHD solves the MHD equations in a comoving coordinate frame. 

g + Vw-o (1) 

— h -V • (pvv + p — BB) = — pv pV$ (2) 

ot a a a 

9E 1 r , d, o 2 B 2 , p 

+ _ V .[v(p + J B)-B(B-v)] = — (p« 2 + — )-^v-V<I> (3) 

Cc ft i *y 

^_l V x(vxB)=--B (4) 

with the equation of state 

E = -pv 2 + + -B 2 (5) 

T 7- 1 2 W 

p = p+^ 2 (6) 

Here, p is the comoving density, p is the comoving gas pressure, v is the proper peculiar velocity, 
B is the comoving magnetic field, E is the total peculiar energy per unit comoving volume, p is the 
total comoving pressure, 7 is the ratio of the specific heats, $ is the proper peculiar gravitational 
potential from both dark-matter and baryons, a = (1 + Zi)/(1 + z) is the expansion factor and t is 
time. 

In this formulation, the comoving quantities that are evolved by the solver are related to the 
proper observable quantities by the following equations: 



Pproper = P * a(t) 3 (7) 

Pproper — P comoving * O (8) 

V proper — ^comoving dx (9) 

^proper = $ - ^adx 2 (10) 

_3 

^proper — B comoving Q> 2 (H) 

It should be noted that the relationship between Im proper and B m m n „,;no in equation[TT]is different 



than that stated in other cosmological MHD codes like iLi et al.l (120081 ). This is due to the additional 
expansion factor that we use in equation [U The proper magnetic field decreases proportional to 
a~ 2 in all formulations of the cosmological MHD equations, but in the formulation we use one half 
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power of a is included as a comoving source term and is due to the redshifting of the photons that 
carry the magnetic field. 

For non-cosmological simulations, the same equations hold, but with with a = l,d = and 
a = 0. This effectively removed each appearance of a from the left hand side, and eliminates the 
terms involving a from the right. For ease of reference, these are: 



dp 
dt 



+ V • (pv) = 



9pv 



+ V • (pvv + p- BB) 



dE 
~dt 



+ V • [v(p + E) - B(B • v)] 



-pV<& 
v • V$ 



dB 

~dt 



V x (v x B) 







(12) 
(13) 
(14) 
(15) 



with the same equation of state, equations [5] and [6j Here, p is the density, p is the gas pressure, 
v is the velocity, B is the magnetic field, E is the total energy per unit volume, p is the total gas 
pressure, 7 is the ratio of the specific heats, <3? is the gravitational potential. The mechanism to 
switch between the two systems of equations will be described in section 12.61 

To solve these equations, we operator split eqns CO)-® into four parts: the left hand side of 
equations J]])-©, the left hand side of equation (jlj), the gravitational acceleration (the two terms 
involving V3>), and the expansion terms (the two terms involving a.) These will be discussed in 
sections l2~6l - 12. 71 In section 12.101 we will discuss the dual energy formulation in Enzo for hypersonic 
flows, and in section 12. Ill we will discuss the Adaptive Mesh Refinement algorithm. We first discuss 
the data structures used to carry all this data in section 12.21 

In the following, we will often have cause to separate the purely fluid dynamical quantities 
p,v,E from the magnetic field B. Unless otherwise noted, 'fluid quantities' will refer to the former 
only. 

For ease of reference, we have supplied a schematic summary of the steps involved in appendix 





2.2. Data Structure 



In Enzo, both parallelism and AMR are done in block decomposed manner. Each patch of space, 
called a grid, is treated as a unique computational problem with Dirichlet boundary conditions 
which are stored in a number of Ghost Zones (see section [231 ) The number of ghost zones depends 
on the method used. The pure-hydro methods in Enzo, ZEUS and PPM, use 3 ghost zones. The 
method we describe here uses 5 ghost zones. 
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Grids are arranged in a strictly nested hierarchy, with each grid having a cell width half that 
of its parent (pure hydro Enzo can take any integer refinement, but the interpolation for MHD is 
restricted to factors of 2.) See figure [TJ Each processor keeps a copy of the entire hierarchy, while 
only one of the processors actual stores the data. 

For all physics modules described in this paper, an individual grid cares not for where it sits 
in space or the hierarchy, and communicates with other grids only through boundary condition fills 
(section ESJ and the AMR cycle (section ETTT]) . 

EnzoMHD in its default mode tracks 14 fields, stored at 3 different points of the cell. The 
5 hydrodynamic quantities, p,v,E to tai are stored at the center of the cell, denoted (i,j,k), and 
represent the volume average of the respective quantities. These are the same quantities stored in 
non-MHD Enzo. 

EnzoMHD tracks 2 copies of the magnetic field and the electric field. One copy of the magnetic 
field is stored in the face of the cell perpendicular to that field component, and represents the area 
average of that field component over that face. This is the primary representation of the magnetic 
field. So Bf^ x is stored in the center of the x face, denoted (i — ^,j,k), Bf, y in the y face at 
(i, j — ^, k), and Bj z in the z face at k — i). It is this field that remains divergence free under 
the cell centered divergence operator: 

V ' Bf = ^ B f^+hi,k ~ B f,x,i-lj,k) + 

'h,,,, •>)• (16) 

The magnetic data structures are one element longer in each longitudinal direction, so for an nx x 
ny x nz grid patch, the Bj x structure is {nx + 1) x ny x nz. 

The second representation of the magnetic field is centered with the fluid quantities at the 
center of the cell. This field is used wherever a cell centered magnetic quantity is needed, most 
notably in the hyperbolic solver in section 12.61 It's equal to the first order average of the face 
centered magnetic field: 

B c^},j,k = °- 5 * ( B f,x,i+±,j,k + B f,x,i-±,j,k) 

B Tyh* = 0-5 * + \,k) (17) 

B SLk = °- 5 * ( B f^,j,k + i + B f, z ,i,i,k-\) 

The final data structure used in EnzoMHD is the Electric Field, which is stored along the 
edges of the computational cell. This represents a linear average of the electric field along that line 
element. Each component is centered along the edge its parallel to, so E x lies along the x edge of 
the cell at — 5, k — |), etc. It is longer than the fluid fields by one in each transverse direction, 
so E x would be nx x (ny + 1) x (nz + 1). 
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Each grid also stores one copy of each of the above mentioned fields for use in assigning ghost 
zones to subgrids. This is described further in 12,51 A temporary field for fluxes is also stored, which 
exists only while the hyperbolic terms are being updated. This data structure is also stored on the 
faces of the zone. There are three fluxes for all 7 MHD quantities. 

For other configurations of EnzoMHD, more or fewer fields may be used. In purely isothermal 
mode (which is at present an option only in EnzoMHD, not in Enzo) the total energy field is not 
tracked, and the isothermal sound speed is taken as a global scalar quantity. This reduces the 
number of fields tracked everywhere the total energy shows up. With dual energy formalism on (see 
section 12 . 10|) an additional field corresponding to either gas energy or entropy is stored, giving an 
additional field where needed. Future work will include multi-species chemistry and more complex 
cooling, which will include additional fields for each species. 

2.3. Consistency 

In several places throughout the flow of Enzo, there may be more than one data structure using 
and writing to a given variable at a given point in space. Ghost zones and face centered fields 
(fluxes and magnetic fields) are examples of this. In EnzoMHD, it is imperative that all data at 
a given point is identical, regardless of the data structure describing it. This may seem like an 
unnecessary comment, but it isn't; in pure hydro simulations, numerical viscosity will damp out 
small perturbations caused by slight inconsistencies in data description. Thus in practice, especially 
in large, stochastic simulations, errors can go unnoticed. Often these discrepancies are negligible, 
other times not, especially when one is concerned with the conservation of a particular variable, 
like V • B. By construction EnzoMHD preserves V • B to machine precision, but it never forces 
V • B = 0; so if it's not zero at the beginning of a time step, it's not going to be at the end, either. It 
is also worth mentioning that inconsistencies in any quantity will cause inconsistencies in the flow, 
which will in turn cause V • B issues. Thus any improper handling of any fluid quantity will cause 
errors in V • B that will persist and usually grow to catastrophic proportions in a relatively short 
period of time. 

There is a prominent redundancy in the magnetic field, namely the field on the surface of the 
active zones of grids. See figure [21 Care is taken to include enough ghost zones, and frequent 
enough ghost zone exchange between grids, that after a time step, two neighboring grids have 
reached exactly the same answer on the surface between the two grids completely independently. 

2.4. Time Stepping 

Enzo uses hierarchical time stepping to determine it's time step. The minimum of 4 different 
criteria is taken for each level, which will be described in sections [2.4.1l - 12~4.41 Timesteps are taken 
in order of coarsest to finest, in a 'W cycle. See figure [3l Given 3 levels, level takes the first step 



- 7- 



of At. Then level 1 takes a single step of At/2. Then level 2 takes one step of At/4. Then, given 
that there are only three levels, it takes another timestep so it is temporally in line with the level 
above. The last three steps repeat: level 1 then takes its second and final step of At/2 so it is now 
at the same time as level 0, followed by two steps on level 2. 

In principle, if a given level has a cell size Ax and the next level of refinement has cell size 
where r is the refinement factor, the more refined grid will have, in principle, time step size —r. In 
Enzo, the step size is chosen for each level and each subgrid time step. In practice, owing to more 
finely resolved structures having slightly higher fast shock speeds, fine grids may in fact take more 
than r time steps for each parent grid step. In some rare cases, such as cosmological expansion 
limiting, a finer grid may take less than r steps. 



2.4-1- Time Stepping: Hydro 

For the hydrodynamics, the harmonic mean of the 3 Courant conditions is used. This was 
demonstrated t o be th e most robust time stepping criterion possible for multi dimensional flows by 



Godunov et al.l (1196 lh . 

Athydro 



l/t x + 1/ty + l/t z 



Ax 

t x =min{ ) (18) 

t v =min( ) 

c f,y 

,Az, 
t z =min{ ) 

C f,z 

where the min is taken over the zones on a level, and Cf jX ,Cf jV and cj jZ are the fast MHD shock 
speeds along each axis: 



2 ^ p 
and similar definition for the other two. 



4. = \ ( a2 + — + J^ 2 + — ) 2 - Aa2B "/p ) ( 19 ) 



2.4-S. Time Stepping: Gravitational Acceleration 

The time step is also restricted to be less than the time it takes for the gravitational acceleration 
alone to move a parcel of fluid half of one zone. 

1 Ax 

At acce i = min(-sqrt ) (20) 

where i = x,y,z and the min is taken of the zones on a level. 
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2.4-3. Time Stepping: Cosmological Expansion 

An additional restriction comes from the cosmological expansion, requiring the timestep to be 
less than the cosmological expansion timescale, 

^t-expansion — (^1) 

where r\ is typically 0.01. 



2.4-4- Time Stepping: Particle Motion 

The fourth timestep criterion is based on restricting particle displacement in a single timestep 
to be smaller than a single zone: 

n V\ x 

particles = min( ) (22) 

^i,p 

where min is over velocity component i and particle p. 



2.5. Boundary Conditions and Ghost Zones 

Ghost Zones are filled in one of three means. 

1. Copying. The dominant mechanism for filling ghost zones copying from active zones that 
occupy the same physical space. This also takes into account periodic boundary conditions. 
For EnzoMHD, face centered fields are copied from the faces of all cells, including those that 
border on active cells. This is somewhat redundant for reasons described in 12.31 

2. External Root grids that lie along the domain wall filled with the external boundary routine. 
If the external boundary condition is not periodic, the grids zones are filled by a predeter- 
mined algorithm; for instance, outflow boundary conditions set ghost zones to be equal to the 
outermost active zone, akin to a Neumann condition of zero slope. These involve outflow, 
reflecting, and a completely general 'inflow'. Note that this is called only on the root grid, 
and not on subgrids that happen to lie on the edge. This can cause spurious waves at reflect- 
ing or outflow boundaries with AMR. Also note for EnzoMHD, the only external boundary 
conditions that have been tested are periodic and outflow. 

3. Interpolation The third mechanism is used on refined grids whose ghost zones do not occupy 
the active space of another grid; these grids have their ghost zones filled by interpolation from 
the parent grid. Since Enzo uses hierarchical time stepping, subgrid steps that begin in the 
middle of a parent grid step fill their ghost zones from a linear interpolation of the parent grid 
time steps at t n and t n+1 . 
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2.6. Left Hand Side: Hyperbolic terms 



With the exception of the l/o term that appears in front of each V- operator, the left hand side 
of equations [D0] are the familiar Ideal MHD equations. A form of equations ((TJ) - (SJ) more relevant 
for this treatment is the following: 

(23) 



dV dF _ 

dt dx 



where 



V 



( P 

PVx 

P Vy 

PV Z 

By 

B z 



(24) 
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B X Vy 

- B, r v y 



\{E + p + B 2 /2) 



p=(E 



-E z 

Ey 

B x (B-v)J 

1)) 



(25) 



(26) 



These form a hyperbolic system of equations, which have been studied extensively in the literature. 
To take advantage of the work already done on this type of system of equations for our cosmological 
algorithm, we first multiply the cell width dx by the expansion factor a. This allows us to use any 
non-cosmological solver for cosmological applications. Upon completion of the solver, dx is divided 
by a to restore dx to the original comoving value. 

Equation [23] is solved by first re- writing it in conservation form, that is taking suitable integrals 
in time and space. The resulting update is, in one dimension, 



At 



p n+ 2 ) 
x,i— | ,j,k' 



(27) 



where V represents the spatial average of the conserved quantities, and F represents an space and 
time average of the flux, centered in time at t = t + At/2. V is the quantity we store in the cells, 
and F comes from the hyperbolic solver. 



The solver we use to solve the hyperbolic equations is that of iLi et al.l (|2008l ). which is comes 
in three parts: spatial reconstruction, time centering, and the solution of the Riemann problem. 
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Spatial reconstruction is done using piecewise linear monotonized slopes on the primitive vari- 
ables (p, v,p, B). Time centerin g of the interface states by A t/2 is performed usin g either the 
MUSCL-Hancock \hi et al.1 [2OO8I ) or Piecewise Linear Method |Colella fc Glazlll98d) inte gration. 
The Rie mann problem is then sol ved using either the HLLC Riema nn solver of IlJ ( 2005 ). HLLD 
solver of iMiyoshi fc Kusand (|2005l ) , or the isothermal HLLD solver of lMignond (|2007l ) . These fluxes 
are computed for the conserved, cell centered variables (p, pv, E, B c ). These fluxes are then differ- 
enced to obtain the update values of the fluid quantities only. The fluxes for the magnetic field are 
stored for use in the Constrained Transport algorithm, discussed in section l2"77l This is done in one 
dimension on successive sweeps along the x, y, and z direc tions. To reduc e operator splitting error, 
the order of the sweeps is permuted. For more details, see Li et al. (|2008l ). 



In isothermal mode, the same method is used, but the energy terms in V and F are removed, 
and only the isothermal HLLD can be used. 



2.7. Constrained Transport and the Divergence of B 



One of the biggest challenges for an MHD code is to ma intain the divergence free constraint 
on the magnetic field (V • B = 0). Brackbill k, Barnes (jl98ol ) found that non-zero divergence can 
grow exponentially during the computation and cause the Lorentz force to be non-orthogonal to the 
magnetic field. There are three major ways to assure the divergence rema i ns ze ro. The first is a 
divergence-cleaning (or Hodge Projection) approach by lBrackbill fc Barnes! (119801) . which so lves an 
extra Poisson's equation to recover V • B = at each time step. But iBalsara fc Kiml (120041 ) found 
that non-locality of the Poisson solver introduces substantial spurious small scale structures in the 
solution. Additionally, solving Poisson's equation on an AMR mesh is computationally expensive. 
Th e second method involves extending the MHD equations to include a divergence wave, as done 
by lPowell et al.1 (|1999l ). iDedner et al.1 (|2002l ). which then advects the divergence out of the domain. 
As most of our solutions are done on periodic domains, this is also an undesirable solution. The 
third method, and the one we have employed in Enzo, is the constrained transport (CT) method of 
Evans fc Hawleyl t) 19881 ) . This method centers the magnetic field on the faces of the computational 
cells and the electric field on the edges. Once the electric field is computed (more on this later) it's 
curl is taken to update the magnetic field. This ensures V • B = for all time, provided it's true 
initially. 



D 



n+l 

f,x,i-\,j,k 



B n . 

X,l- 



2 >J ~ 9 > 



E_ 



Az 



,)+ 
)) 



(28) 



Plugging equation IB12I into the divergence operator [16] to find V • B^ +1 , one finds all terms are 
eliminated except the initial divergence V • Bj. 

The CT algorithm of Evans fc Hawley (|l988l ) was extended to work with finite volume methods 
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by iBalsara h Spicerl (|1999l ). This method uses the fact that the MHD Flux has the electromotive 
force as two of its components (see the 5 th and 6 th components of eqn. I25j) , so using these components 
then incorporates all the higher order and shock capturing properties of the Godunov solver into the 
evolution of the electric field. These components, which are centered at the face the computational 
cell, are then averaged to obtain an electric field at the edges of the cell. This was the first CT 
method applied to Enzo, so unless otherwise n oted, the simulations pre sented here were done with 
this method. The reader is encouraged to read IBalsara &: Spicerl (|1999i ) for the full details. 



Gardiner Stond (|2005l ) extended this idea to in clude higher order spati al averaging, which 



eliminates a number of numerical artifacts present in IBalsara Spicerl (|1999l ) and increases the 
accuracy of the method. This method uses the fluxes from the Riemann solver, plus additional 
information from the data in the cell to construct a linear interpolation from the cell face to the cell 
edge. The reader is encouraged to see that paper for the details. 

After the curl is taken and the face centered field Bf is updated, it is then averaged to obtain 
B c , via equation 1171 



2.8. Right Hand Side: Gravitational Acceleration 

In cosmological simulations, Enzo tracks the proper peculiar gravitational potential. 

V 2 $ = —(p b + Pd-po) (29) 
a 

where pt, and p<i are baryonic and dark matter comoving density respectively, and po is the comoving 
background density. For non-cosmological simulations, the dark matter and background density are 
ignored. 

The gravitational potential <3? is solved in Enzo using a combination of methods. First, the 
root grid potential (which covers the entire computational domain) is solved for using a fast Fourier 
transform. Then the subgrids (which hopefully do not cover the computational domain) are solved 
using a multigrid relaxation technique. This resulting potential $ is then differenced to obtain the 
acceleration g = V<3?. Specifically, 

g« = -(*»+! - *<-i) (30) 

As mentioned before, the fluxes are computed at the half time point t + l/2At. In order to 
keep the velocity and consistent with this time centering, they are first advanced by a half time 
step: 

At , , 

v = v + — g (31) 

After the fluxes are differenced to obtain the new state , these states are then updated 
with the accelerations. For the velocity update, a density field centered in time is used. We follow 
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the same formulation used by IColella &z Woodward! (|1984l ) 

v x n ^ - 

1 



P 



,n+l 



E 



n+l 



E 



/n+l 



n+l/ , /n+l\2 



(V. 



(32) 
(33) 



2.9. Right Hand Side: Expansion Source Terms 

The cosmological expansion source terms are treated in much the same way as the gravitational 
source terms. First, a half time step is added to the values before the flux is computed. 



v 'n =v n _ i^f* n (34) 

2 a 

p'n =p n _ I A^3( 7 - l)p n (35) 

B? =B n c - l^B n c (36) 



The quantities v' n , p' n and B' n are then used in the rest of the solver described in section l2~6l 
After the fluxes are differenced, the source terms are then added to the fluid quantities in full. This 
is done in a semi-implicit manner, by averaging the quantities to be updated in time. For instance, 
the expansion contribution to the magnetic field is 

dB a 



which is discretized 



* 2a B (37 » 



gn+1 _|_ gn+l\ 

B n + p l - B n+1 = _iL(-S2— -) (38) 



and solving for B n+1 we have 

exp 



2a 



a 

>n+l _ (j- ~ x ) pn+l 



(39) 



exp (1 + x) v ' 

Pressure and velocity are updated in a similar manner. See appendix [O for the full update. 



2.10. Dual Energy Formalism 



Hypersonic flows are quite common in cosmological simulations. Due to the extremely large 
gravitational forces, the ratio of kinetic energy Ek inetic to gas internal energy E interna i can be as 
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high as 10 s . This leads to problems when computing the internal energy in this type of flow, as 
the universe does math with infinite accuracy, but computers do not. Higher order Godunov code 
typically track only the total energy (equation [5]). Thus finding the internal energy from the total 
energy tracked by the software, 



^internal — ^total Ej* 



inetic 



E 



magnetic 



involves the small difference of two (or three) large numbers, which causes problems when the small 
number (E interna i) is near the roundoff noise of the original numbers (E tota i and E kinetic + E magnetic ) . 

To overcome this, we have implemented two algorithms th at solve an addit ional equation to 
track the small numbers; the modified entro py equation given in lRyu et al.l (|1993l ) and the internal 
energy equation given in 



Bryan et al.l t| 19951 ) . These two equations are: 

3( 7 - Da 



dS 1„ ,„ s 



-s 



3(7-1)6 



pe + -V 

a 



(41) 
(42) 



where S = p/p 1-1 is the comoving modified entropy and e is the internal energy The modified 
entropy equation is valid only outside the shocks where the entropy is conserved. Use of either (not 
both) of these equations is at the discretion of the simulator. 

Through the course of the simulation, the ratio of internal energy to total energy is monitored. 
When this rat io is less than some preset value r], one of the modified equations is used. As in 
Li et al.l (|2008l ). we use rj = 0.008. They note that reducing this parameter will cause a decrease 
in the volume filled by low temperature gas, as most of the gas affected by the switch is cold, high 
velocity gas. T he optimal choice for this parameter is still an open question for the general situation. 
Li et al.l (120081 ) compared this two approaches and found almost identical results. 



2.11. Adaptive Mesh Refinement 



Structured AMR, initially devised by Berger fc Colella (jl989l ). is a technique for increasing res 



olution of a simulation in parts of a simulation that require higher resolution for increased accuracy 
or suppression of numerical artifacts, while conserving memory and CPU cycles in areas that don't. 
Refinement criteria will not be described here, as they vary from simulation to simulation. AMR 
has four basic necessary parts: 

1. Patch Solver This is the algorithm that actually solves the finite volume PDEs in question, as 
described by sections l2~6l - 12.101 The approximations used for the patch solver are conservative 
in a finite volume sense, and the rest of the choices are made to preserve that conservation. 



2. Refinement Operator This is the routine that creates fine resolution elements from coarse 
ones. In Enzo, we use conservative, volume weighted interpolation for the fluid quantities 
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p,E,v. For the magnetic fields, we use the method described by iBalsaral (|200ll ). with some 
slight modifications in implementation. This method constructs a quadratic divergence free 
polynomial, and area-weighted averages are used for the fine grid quantities. This is described 
in more detail in appendix lAl 

3. Projection Operator This is the routine that projects the fine grid data back to the parent 
coarse grid. For Enzo, the parent grid is simply replaced by a volume-weighted average of the 
fine cells. For the face centered magnetic field, this is an area weighted average, though in 
practice we don't explicitly average the magnetic field, as discussed in below and in appendix 

4. Correction Operator Once the projection operator replaces the solution on the coarse grids, 
the evolution on the coarse grids is no longer consistent with the underlying equations in the 
manner they were discretized. That is to say, the total change of any conserved quantity inside 
the region is no longer equal to the flux across its surface. For the Enzo hydro fields, this is 
corrected with the flux correction mechanism. More details on this and the modifications in 
EnzoMHD see appendix iBl 



EnzoMHD does all of these steps for the fluid quantities, but for the magnetic field it slightly 
alters this procedure. In order to overcome a shortcoming in the original data structures used in 
Enzo, we combined the projection and correction operations for the magnetic fields in one step. 
The net effect of the correction operator is to ensure that all zones are updated by finest resolution 
fluxes available, even if they were updated by coarse data initially. For the magnetic field update, 
we don't project the actual magnetic field that is of interest, but rather the electric field (effectively 
the 'flux' for JB/), then take the curl of the newly projected electric field. Thus the coarse magnetic 
data co-located with the fine grids get updated with the fine data, and the bounding zones don't 
need correction at all. 

More detail on this process can be found in appendices [A] and [B] 
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Fig. 1. — A schematic of a parallel AMR hierarchy on two processors (l eft) and a grid patch with 
ghost zones (right). Image courtesy James Bordner, initially appeard in (jNorman et alJl2007l ) 
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Fig. 2. — Data redundancy of the face centered magnetic fields: the face centered field denoted by 
the stars are updated by both grid 1 and grid 2. Enough ghost zones are exchanged to ensure that 
the entire stencil for the update of these fields is the same in both data structures. 
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Fig. 3. — A depiction of the timestep strategy in Enzo 
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3. Numerical Experiments 

EnzoMHD has many configurations available. Here, we test some of the possible configurations, 
to indicate the quality of solution possible with EnzoMHD. 



3.1. MHD Tests without AMR 

We first test our code in unigrid (fixed resolution) mo de, in order to ensure consistency of the 
patch solver with the algorithm described in Li et al. (hood ). We do two one dimensional cosmology 



tests (Caustics and Zel'dovich Pancake), two one dimensional non-cosmological tests (Brio and Wu 
and the Kim Isothermal), one 2d non-cosmological test (Orszag Tang) and one 3d cosmological test. 



3.1.1. Brio and Wu shock tube 



The shock tube defined by iBrio h Wul (|1988l ) is a standard test of any MHD solver, as it 



displays a number of the important MHD waves, including a compound wave. Compound waves 
are not a property of pure hydrodynamics, because the system is convex. However, due do the more 
complex nature of the MHD equations, certain initial conditions can cause flows in which at one 
point the shock speed in a given family is higher than the wave speed for that family, causing a 
shock, but lower in the post shock region, causing a rarefaction immediately following the shock. 

This can be seen in figured! The problem was run with 800 zones to a time t = 0.2, using the 
HLLD solver in Enzo. This shock tube shows, from left to right, a fast rarefaction, slow compound 
(shock+rarefaction), contact, slow shock, and fast rarefaction. It can be seen that this solver 
captures this shock tube problem quite well. 



3.1.2. Isothermal Tests 

One of the primary application areas of EnzoMHD will be in simulating turbulence and star 
formation in cold molecular clouds. Due to the fast cooling time of these environments, an isothermal 
eq uation of state is a g ood approximation a large portion of these processes. In simulations done 



by iKritsuk et al.l ([2007]) using Enzo and other works by the same authors an isothermal equation of 



state is approximated by using an adiabatic solver and setting 7 = 1.001. 

To tes t if th is approximation is appropriate for this code, we ran the isothermal shock tube of 



Kim et al.l t| 19991 ) . One can see from figure that this approach works well, as shock jumps and 
positions are all correct, and features are reasonably sharp. This test was run with 256 zones to a 
time of 0.1. 



However, in turbulent simulations with gravitational collapse, the measured value of the sound 
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Fig. 4. — The shock tube of lBrio h Wul (119881 ) . showing from left to right a fast rarefaction, slow 
compound (shock+rarefaction), contact, slow shock, and fast rarefaction. T=0.08, and 800 zones 
were used. 
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speed, yjp/p, is initially uniform, but after a few hundred timesteps can vary by as much as 1000, 
which is far from isothermal. It is believed that the difference between this code and what has been 
done in the past with Enzo stems from the Riemann solver. The HLL family of Riemann solvers 
assumes a particular wave structure in computing the interface flux. This wave structure, for HLLC 
and HLLD, contains a contact discontinuity which is not present in the isothermal Riemann fan, 
and does not reduce appropri a tely i n the 7 — ► 1 limit. To combat this, we installed the Isothermal 
variant of HLLD by iMignond (|2007l ). The results of this code on the Kim test are nearly identical 
to that in figure [5] and not reproduced here. The problem seen are, of course, eliminated as the 
sound speed is set as an input parameter. 



3.1.3. One- dimension MHD Caustics 



This test is taken from lLi et al.l (|2008l ). which initially derived from a pure hydro version from 
Ryu et al.l (119931 ). This problem is used to test the ability of the code to capture shocks and to 
deal with hypersonic flows. Initially, v x = —^sin(2irx) 1 p = 1 and p = 10~ 10 . Caustics are formed 
because of the compression by the velocity field. The Mach number of the initial peak velocity is 
1.2 x 10 4 . The pressure can easily become negative for such high Mach number flow. 



We performed the test with same magnetic field settings as in lLi et al.l (120081 ). The magnetic 
field in the x and z directions are always zero while B y = 0,0.001,0.02 and 0.05. The calculation 
was done with 1024 cells and the results at t = 3 are shown in figure [6l Our results match the 



results from CosmoMHD (ILi et al .1120081 ) quite well, as expected 



3.I.4. The Zel'Dovich Pancake 



The Zel'Dovich pancake is a popular test pro blem for codes t hat include gravity in comoving 
coordinates. The problem setups are taken from Li et al. (j2008l ). This takes place in a purely 
baryonic universe with Q = 1 and h = ^. The initial scale factor ai = 1 corresponds to Zi = 20. 
The initial velocity field is sinusoidal with the peak value 0.65/(1 + Zj), and v = at the center 
of the box. The initial comoving box size is 64/i _1 Mpc. The shocks forms at z = 1. The initial 
baryonic density and pressure are uniform with p = 1 and p = 6.2 x 10 -8 . The tests were run with 
1024 cells, both wi th and withou t magnetic fields. Our results are almost identical to the results 
from CosmoMHD (|Li et al.ll2008l ). as expected. Results can be seen in figure [7l 



3. 1 . 5. Orszag- Tang 



The Orszag- Tang Vortex was originally developed by lOrszag fc Tang) (Il979l ) to demonstrate 
that small scale structure can be generated by the nonlinearities in the MHD equations. It initially 
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Fig. 6. — 1-D MHD caustics at t = 3. Density, gas pressure, total pressure and B y are plotted. For 
the small field runs, almost no change can be seen, while larger field runs decrease the peak of the 
density considerably due to the increased pressure. 
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Fig. 7. — The Zel'Dovich Pancake problem with various values of the magnetic field, at t = 0. 
Increasing the magnetic field strength increases the central magneti c pressure, redu cing the density 
and changing the overal solution structure. Results match those of iLi et alJ (j2008l ). 
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starts with a single large scale rotating velocity structure and two circular magnetic structures. 
From these simple large scale initial conditions, substantial small scale structure is formed. It now 
serves as a standard test problem to demonstrate the accuracy and diffusivity of MHD codes. 

The initial conditions are on a 2 dimensional periodic box, 256 zones on a side, v = Vq( — sin{2iry)x+ 
sin(27rx)y,~B = Bq{— sin{2%y)x + sin(4irx)y), Vq = 1,Bq = 1/a/47t, Po = 25/(367r),p = 5/(127r), 
and 7 = 5/3 which gives a peak Mach number of 1 and peak [5 = Pq/{Bq/2) = 10/3. Figure [8] 
shows the density at t = 0.48, from which one can see that the solution agrees with other solutions 
to the problem in the literature. 



3.1.6. 3D Adiabatic Universe with MHD 



We have also performed the 3D adiabatic CDM Universe test described by Li et al. (I2OO8I ) both 
with and without magnetic fields. We also compared the non-magnetized results with the results 
run using the PPM solver (|Colella fc Woodward! 11984 ) . Adiabatic evolution of a purely baryonic 
Universe was computed with an initial CDM power spectrum with the following parameters: £2 = 
Qb = 1, h = 0.5, n = 1 an d <7s = 1 in a c ompu tational volume with side length L = 64/i _1 Mpc. 
The transfer function from bardeen et al.1 (jl986l ) was used to calculate the power spectrum of the 
initial density fluctuations. Evolution was done from z = 30 to z = 0. We used 256 3 cells for each 
simulation. The comparisons are made at the final epoch, z = 0. Though this test is identical to 



that of lLi et al.l (120081 ). our results can't compared with theirs directly since different random seeds 
were used for the realization of the initial density and velocity. 

Figure M shows a comparison of the mass-weighted temperature distribution, figure [10] is a 
comparison of the volume-weighted density distribution. The discrepancies between PPM and 
MHD solvers are small, indicating the two codes perform roughly the same. The nature of the 
differences is expected, since PPM solver has third order accuracy while the MHD solver has second 
order accuracy and larger numerical diffusion. This allows PPM to capture shocks in fewer zones, 
which causes the dense shocked gas to not only have a smaller volume fraction, but also be hotter 
and slightly less dense than in the MHD solver. 

We have also done a similar run with the same initial conditions to the above, but with an initial 
magnetic field, B x = B z = 0, B y = 2.5 x 10~ 9 Gauss, which is 4.32 x 10~ 7 in code units. Figure 
[TT] shows the scaled divergence of the magnetic fields, averaged over the entire box, as a function of 
redshift. The scaled divergence is < |/iV-B/|B|| >, where h = 1/256 is the spatial scale, and \B\ is 
the local maximum magnetic field strength, is the most relevant measure of the potential numerical 
effects of divergence. The divergence of the magnetic fields is close to the round-off error. 



Fig. 8. — Density from the Orszag-Tang vortex, at t=0.48. Initial conditions are uniform density, 
with a single rotating velocity structure and two circular magnetic structures. This generates 
significant small scale structure, which has been used to compare effective resolution of different 
MHD schemes. 
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Fig. 9. — Comparison of mass- weighted temperature histogram at z = for the 3D purely baryonic 
adiabatic Universe simulation. The solid line is from the MHD code and the dashed line is from 
Enzo-PPM. 
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Fig. 10. — Comparison of volume-weighted density histogram at z = for the 3D purely baryonic 
adiabatic Universe simulation. The solid line is from the MHD code and the dashed line is from 
Enzo-PPM. 
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Fig. 11. — The scaled divergence < \hV • B/|B|| > of the magnetic fields for the 3-D simulations 
of a purely baryonic adiabatic Universe. Here h = 1/256 is the scale length and ||5|| is the local 
magnetic field norm, and the average is over the entire volume. Scaled divergence is a more relevant 
measure of numerical effects of divergence than the strict divergence. As desired, the divergence is 
near the machine round off noise, the theoretical limit. 
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3.2. 



MHD Tests with AMR 



To test the Adaptive Mesh Refinement, we ran a sample of the tests from the previous section 
with AMR, to ensure no spurious artifacts are introduced by the AMR. These are the Adiabatic 
Expansion test in section 13.2.11 and the one dimensional caustic and pancake tests (sections 13.2.21 



This test is taken from iBryan et al.l (119951 ). This test uses a completely homogenous universe 
with initial T\ = 200K and v% = WOkm/s in the x-direction at an initial redshift of z% = 20. In 
the code units, the initial density is 1.0 and initial velocity is 2.78 x 10~ 3 and the initial pressure 
is 1.24 x 10 -9 . Additionally we have a uniform magnetic field B x = B y = B z = 1 x 10 -4 in code 
units, which is 2.66 x 10~ 7 G in cgs. 

The simulation used a 16 3 root grids with 2 levels of refinement in the center region and ran 
to z = 0. 

The expansion terms in eqns - (SJ) operate like drag terms, so that in the absence of a 
source, the velocity decreases as v = Via -1 , the temperature as T = Tia~ 2 and the magnetic field 
should decrease as a~ x / 2 . 

The temperature at z = is 0.453406K, 0.024% below the analytic result of 0.453515K. 
The velocity at z = is 4.76176fcm/s, compared to the analytic result 4.7619fcm/s, a 0.0029% 
discrepancy. The final magnetic field strength is 6.03 x 10 _10 G (2.18 x 10~ 5 in the code units), a 
difference of 0.0006% with respect to the analytic solution. Figure [12] shows the B y as a function 
of redshift, the solid line shows the theoretical value. 



We also ran the the Id MHD Caustic test with AMR, using 256 root grid zones with 2 levels 
of refinement, again by a factor of 2, giving an effective resolution is 1024 cells. Figure [T3l shows 
comparisons of density and gas pressure of non-AMR and AMR runs with different initial magnetic 
field strengths, as described before. Figure E] shows the comparisons of B y for runs with different 
initial values of B y . In both plots, the AMR result is sampled to the finest resolution. The AMR 
runs give almost identical results to the unigrid runs, while the CPU time and memory were greatly 
saved in the AMR runs. 



and Ed- 



3.2.1. Three- dimension MHD Adiabatic Expansion 




3.2.2. One- dimensional MHD Caustics with AMR 



- 29 - 




Fig. 12. — Magnetic field in the y direction in the AMR MHD adiabatic expansion test. The pluses 
show the results of simulation and the solid line is the analytic result. 




Fig. 13. — Comparisons of density and pressure in the MHD Caustic tests, non-AMR vs AMR. The 
left column shows density and the right column shows gas pressure. Initial magnetic field of each 
row from top to bottom is 0, 0.001, 0.02 and 0.05. 
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Fig. 14. — Comparison of B y in the MHD Caustic tests, non-AMR vs AMR. Initial magnetic field 
of each panel from top to bottom is 0.001, 0.02 and 0.05. 
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3.2.3. Zel'Dovich Pancake with AMR 

We also ran the pancake problem with AMR. The problem was set up with the same initial 
conditions as the unigrid run, but with a root grid of 256 root cells and 2 levels of refinement by 2. 
We compared these results having effectively 1024 cells to the results of our previous high resolution 
which actually had 1024 cells. Figure [11] shows comparisons of density and gas pressure between the 
non-AMR and AMR runs, with different initial values for B y . Figure [16] shows the comparisons of 
By with different initial values. Again, the AMR computation got very similar results, while saving 
CPU and memory resources. 



3.2.4- MHD Galaxy Cluster Formations 



Cluster fo r mation (without MHD ) has been studied intensively by researchers using Enz o 
|Normanl 1200.4 brvan fc Normanl Il998l : lLoken et all I2OO2I : iMotl et ali I2OO4I : lHallman et alJbood ). 
It is one of the most important applications of Enzo's high dynamic range. Many cluster simula- 
tions have been run with Enzo with a wide variety of physics (i.e. radiative cooling, star formation, 
etc) and we can compare these results to similar simulations run with MHD. More information about 
Enzo simulated cluster can be found in Simulated Cluster Archive at http:/ /lea. ucsd.edu/data/sca/. 
Here, we present just one simulation to demonstrate the MHD code. 

This simulation uses a Lambda CDM cosmology model with parameters h = 0.7, Q m = 0.3, 
= 0.026, J7a = 0.7, as = 0.928. The survey volume is 256 h~ x Mpc on a side. The simulations 
were computed from a 128 3 root grid with 2 level nested static grids in the center where the cluster 
form. This gives an effective resolution of 512 3 cells (0.5 h~ l Mpc per cell) and dark matter particle 
mass resolution of 1.49 x 10 10 solar masses initially in the central region. Adaptive mesh refinement 
is allowed only in the region where the galaxy cluster forms, with a total of 8 levels of refinement 
beyond the root grid, for a maximum spatial resolution of 7.8125 h~ x k"pc. While the baryons are 
resolved at higher and higher spatial and mass resolution at higher levels, the dark matter particles 
maintain constant mass so as not to add any additional noise. The simulations are evolved from 
z = 30 to z = 0, and all results are shown at the redshift z = 0. We concentrate our study on a 
cluster of M = 1.2 x 1O 15 M . 

In order to isolate the effects of the numerical approximation from the effects of MHD, we first 
run the simulations adiabatically without additional physics and the magnetic field set to zero, and 
compare to a PPM run with identical parameters. In table [H we list the basic parameters for the 
clusters formed in each solver. The viral radius, R v i r is calculated for an over density S -f of 200. 
M v i r , Mdm and M gas are the total mass, mass of the dark matter and mass of the baryon inside 
the virial radius, respectively. T v i r is the average of the temperature of the ICM inside the virial 
radius. Evidently, there is very little difference between the results from the two solvers. 



Figures [l~7lfl~9l show the images of the logarithmic projections of the dark matter density, gas 



33 



■ nmniiiiiiiiiiiiiiiiiiiiiiiiii n iiiiiiii ^ | AMR 



+ | ^ IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIM 




0.42 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 0.42 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 



42 0.44 46 0.48 Oh 0.52 0.54 56 0.58 0.42 0.44 46 48 5 52 54 56 58 




■ iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii n l 




- i iiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiii- 



0.42 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 0.42 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 




0.42 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 0.42 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 



Fig. 15. — Comparisons of density and pressure in non-AMR and AMR runs of the Pancake test. 
The left column shows density and the right column shows gas pressure. Initial magnetic field of 
each row from top to bottom is 0, 1.3e-6G, 2e-5G and le-4G. 
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Fig. 16. — Comparisons of magnetic y component in non-AMR and AMR runs of the Pancake test. 
Initial magnetic field of each panel from top to bottom is 1.3e-6, 2e-5 and le-4G. 
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density, and X-ray weighted temperature, respectively, at z = 0. Both PPM and MHD solvers show 
very similar images in all three quantities, differing only slightly in the small scale details. 

Figure I20ll22l show the radial profiles of dark matter density, gas density, and x-ray weighted 
temperature. The profiles match quite well in all three quantities, with only minor differences. 
There is a slight deviation in the radial profiles of dark matter density near the center of the cluster, 
but this is near the resolution limit of the simulation, so not a trustworthy data point. In the density 
profile, it can be seen that the MHD solver gives a slightly higher average density. The temperature 
agreement is good enough to not worry about. 

We have also run the simulations with non-zero initial magnetic field. A uniform initial magnetic 
field of 9.72753 x 10~ 10 G (1 x 10" 7 i n code units) in the y direction was added to the system at the 



start of simulation at z = 30. Since iDolag et al.l (| 19991 ) has shown that the initial magnetic fields 



structures are not important to the final magnetic fields structures in their MHD SPH simulations, 
no other initial magnetic fields configuration will be used in this paper. Figure [23] shows 4 projections 
of the cluster center: gas density, temperature, magnetic energy, and synthetic Faraday rotation 
measurement RM = 2 T T m' 2 c i Jo n eBds. We can see that the gas density and temperature images 
are almost identical to the MHD run with zero magnetic fields. As expected, the magnetic energy 
is concentrated in the cluster core. The maximum magnetic fields is 1.0630270 x 10~ 8 C The RM 
is about 2-3 radm~ 2 at the cluster core. Figure [24] shows comparison of the radial profiles of the 
simulations with and without initial magnetic fields, while figure [25] depicts the volume weighted 
averaged radial profiles of the magnetic field strength and plasma f3. Since f3 is quite large, these 
small magnetic fields acts as a passive tracer of the plasma and has little effects on dark matter and 
gas dynamics. 

To further test our code, we also ran a simulation with a relatively large initial magnetic fields. 
We also included radiative cooling, star formation, and stellar feedback. The radiative cooling 
models X-ray line and bremsstrahlung emission in a 0.3 solar metallicity plasma. The star formation 
model turns cold gas into collisionless star particles at a rate psf = Vsf max ( T Pb L Td — y, where tjsf 
is the star formation efficiency factor 0.1, and r coo / and Td yn are the local cooling time and free 
fall time, respectively. Stellar feedback returns a fraction of stars' rest energy as thermal energy at 
a rate Tsf = VsnPsfc 2 to the gas. We did two runs, one without initial magnetic fields and the 



Table 1: Cluster Properties 



Parameter 


Hydro PPM 


MHD 


R vir (Mpc) 


2.22946 


2.22674 


M vir (M Q ) 


1.26462e+15 


1.25999e+15 


M dm (M Q ) 


1.09746e+15 


1.09683e+15 


M gas {M & ) 


1.67158e+14 


1.6316e+14 




8.68422e+07 


8.66301e+07 
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Fig. 17. — Logarithmic projected dark matter density at z = 0. The images cover the inner 4 
Mpc/h of cluster centers. The left panel shows the result from the PPM solver and the right panel 
shows the result from the MHD solver. The color bar is in Mq Mpc~ 3 . 
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Fig. 18. — Logarithmic projected gas density at z = 0. The images cover the inner 4 Mpc/h of 
cluster centers. The left panel shows result from PPM solver and the right panel shows result from 
MHD solver. The color bar is in M & Mpc~ 3 . 




Fig. 19. — Logarithmic projected X-ray weighted temperature at z = 0. The images cover the inner 
4 Mpc/h of cluster centers. The left panel shows result from PPM solver and the right panel shows 
result from MHD solver. The unit is Kelvin. 




Fig. 20. — Spherically averaged dark matter density radial profile at z = from MHD solver and 
PPM solver. 
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Fig. 21. — Spherically averaged gas density radial profiles at z = from MHD solver and PPM 
solver. 
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Fig. 22. — Spherically averaged temperature radial profiles at z = from MHD solver and PPM 
solver.. 
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Fig. 23. — Images of gas density (Mq Mpc~ 3 ), temperature (K), magnetic energy density (erg 
cm~ 2 ) and rotation measure (rad m~ 2 ) of the galaxy cluster simulation with an initial magnetic 
field B y = 9.72753 x 10~ 10 G. Projections are of the inner 4Mpc/h of cluster center at z = 0. 
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Fig. 24. — Specially averaged radial profiles of dark matter density, baryon density and temperature 
of MHD simulations with zero and By = 9.72753 x 10~ 10 G initial magnetic fields. 
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Fig. 25. — Spherically averaged radial profiles of magnetic field strength and plasma f3 of MHD 
simulation with B y = 9.72753 x 10 10 G initial magnetic fields. 
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other is with a large initial magnetic fields of B y = 1.0 x 10~ 4 in code units (9.72753 x 10 -7 Gauss.) 
Figure [26] shows the radial profiles of gas density and temperature of both runs and the magnetic 
field strength and the plasma f3 of the run with magnetic fields. 

The magnetic field s reached 20 fi G in the core region, a few times larger than the observations 
(jCarilli fc Taylor! |2002j) . In the center where (3 reaches a minimum, t he kin etic energy is a few 
percent of the thermal energy, as expected from Iapichino fc Niemeyer (|2008l ). The magnetic field 
has become dynamically important in the cluster center. The effect is not significant in the density, 
as seen in the upper right plot in figure [26j but definitely noticable in the temperature field, as 
some of the thermal pressure that was balancing the collapse is replaced by magnetic pressure. In 
this way, magnetic fields may help to cool cluster cores, giving a better match to observations. 
Detailed analysis of the magnetic field structure and their influence on the cluster will be presented 
in forthcoming paper. 
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Fig. 26. — Radial profiles of MHD simulations with zero and B y = 9.72753 x 10 1 G initial magnetic 
fields with radiative cooling, star formation and stellar feedback. 
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4. Conclusion 



In this work, we have presented the implementation of MHD in the AMR cosmology code Enzo 
in order to serve as a single complete reference document for future simulations done with EnzoMHD, 
and a reference for future users of the code. EnzoMHD is capable of multi-resolution cosmologi- 
cal and non-cosmological astrophysical simulations using ideal MHD. Enzo uses block structured 
AMR, which solves they hydrodynamic (and now magnetohydrodynamic) PDEs on fixed resolution 
patches, and communicates the finest resolution information between coarse and fine patches in way 
that is conservative in the volume- averaged quantities. This entails 4 basic components: the PDE 
patch solver, creation of fine grids (interpolation), communication of fine data back to coarse data 
(projection) and correction of the interface between coarse and fine grids (flux correction). MHD 
has the additional constraint that the divergence of the magnetic field, V • B, must be zero at all 
times, which requires additional machinery to advance the PDEs (Constrained Transport) and some 
modifications to the projection and flux correction steps. In addition to multi-resolution hydrody- 
namics, EnzoMHD includes the effects of gravitational acceleration and cosmological expansion, 
and a modification to the base PDE solver to account for flows with large disparity between kinetic 
and thermal e nergies (dual energy formalism). In EnzoMHD, we used we use the PDE solver of 
Li et al. fj200sl ) to solve the ideal MHD equations (section [276]) for the patch solver, which is second 



order accura te in both tim e and space. We use a slightly modified version of the AMR algorithm 
procedure of iBalsaral (|200ll ) to create interpolate fine grids and project the more accurate fine grid 
da ta to the cheaper c oarse grid data (section 12. Ill and ap pendix lA)) . We have used the CT methods 
of balsara k Spicerl (llQQflh and bardiner k Stone! (120051 ) to advance the induction equation while 
maintaining the constraint V B = Ofsection 12.71 We have operator split the gravi tational (|2.8p and 
cos mological expansion (|2.9p terms; and included the dual energy techniques of iRyu et al.l (119931 ) 
and iBryan et al.l (119951 ) . 



In section [31 we present the results of a broad array of tests to demonstrate the accuracy of the 
chosen methods. These include the shock tube of Brio and Wu 13.1.11 the isothermal shock of Kim 
13.1.21 on dimensional MHD Caustics l3.1.3l the famous Zel'Dovich Pancake 13.1.41 the Vortex problem 
of Orzag-Tang 13.1.51 an adiabatic expanding universe 13.1.61 Some of these were additionally run 
with AMR, and the results compared to the unigrid case. The results of these overall agree with 
both what's been present in the literature before and comparisons with our existing PPM solver. As 
an example of the capability and application area of this code, we present some preliminary results 
from a calculation of galaxy cluster formation with magnetic fields in section 13.2.41 

Currently underway are simulations involving protostellar core formation, MHD Turbulence, 
and galaxy cluster formation and evolution with magnetic fields. Work has begun to include cosmic 
ray acceleration, sink particles for star formation, and ambipolar diffusion into the code. 

This work has been supported in part by NSF grants AST-0708960 AST-0808184, AST-0807768 
and by NASA grant NNX08AH26G. Additional support was supported by IGPP at Los Alamos 
National Laboratory. Simulations described in this paper were performed at the San Diego Super 
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Computing Center with computing time provided by NRAC allocation MCA98N0202 and LANL 
Institutional HPC clusters 



5. Appendix 

A. AMR MHD Reconstruction 

A.l. MHD Reconstruction 

For completeness, we will briefly outline the AMR re construc t ion us ed in EnzoMHD. The reader 
is encouraged to see the details in the original paper by iBalsaral (|200ll ). 



In this appendix, we have dropped the subscript / from the face centered fields, as the face 
centered field is the only one in question. 

Balsara's reconstruction method for the magnetic field is a 3 dimensional, quadratic recon- 
struction of all 3 vector fields simultaneously. If we let b be the polynomial fit to the discrete face 
centered field field -B, the general reconstruction is 



b x (x, y, z) = a + a x x + a y y + a z z + a xx x 2 + a xy xy + a xz xz 
b y (x, y, z) = b + b x x + b y y + b z z + b xy xy + b yy y 2 + b yz yz 
b z (x, y, z) = c + c x x + c y y + c z z + c xz xz + c yz yz + c zz z 2 

The coefficients are found by the following constraints: 



(Al) 
(A2) 
(A3) 



1. The analytic reconstruction should be divergence free. 

2. At the faces of the parent cell, the reconstruction should reduce to a bilinear reconstruction, 
where the slopes are monotonized with the minmod slope limiter. For instance, 



b x{x = —,v) 



B 



x,i+Lj,k 



Ay 



-y + 



Az 



where 



minmod{B x i 



minmod(x, y) 



x, 

y, 



J' B x,i+±,j 

\x\ < \y\ and xy > 
\y\ < \x\ and xy > 



0, xy < 



(A4) 



(A5) 



(A6) 
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The minmod slope is used in order to minimize oscillations. Area weighted averages over these 
polynomials are then used to assign the fine grid values. 

Often, a fine grid patch will encroach on unrefined territory. This results in the refinement 
of coarse zones that a.) share a face with fine grids but b.) don't have corresponding fine grids 
of their own. Balsara refers to this as "Prolongation" of the fine grid. To avoid generating any 
divergence at the boundary of the face, the interpolation polynomials need to match the old fine 
data. The interpolation equations above (eqns [All - IA3|) do not have enough degrees of freedom to 
accommodate that many data points. In this case, Balsara describes a new polynomial that DOES 
have enough degrees of freedom, by adding 3 rd order cross terms to equations I All - lA3t 



b x (x, y, z) =a + a x x + a y y + a z z + a xx x 2 + a xy xy + a xz xz 

~i~ dyzl/Z ~i~ O'xyzXyZ + & XX z% Z ~\~ 0, xxy X y (A7) 

b y (x, y, z) =b + b x x + b y y + b z z + b xy xy + b yy y 2 + b yz yz 

+ b xz xz + b yyz y 2 z + b xyz xyz + b xyy xy 2 (A8) 

b z (x, y, z) =c + c x x + c y y + c z z + c xz xz + c yz yz + c zz z 2 

+ c xy xy + c yzz yz 2 + +c xzz xz 2 + c xyz xyz (A9) 

The yet undetermined coefficients are found by matching the polynomial to a bilinear fit on 
the face: 

Ax . A y £ i+ i A Z B i+ i A yz B xi+ i 

b ^ X = —^ = + ~^y-^ y + -^z~^ Z + ^y~Kr yZV (A10) 

and now the finite differences are taken from the finest grid: 

A yz B X)i+ i = 4((5 a . i i+ i J+ i M i ~ B X)i+ i tj _i M i)- 

(B Xji+ y + i >k _i - B X j + y_i tk _i)) (All) 

A v B x,i+± = ((-B Xji+ i J+ i jfc+ i - 5 a . >i+ i J ._i )fc+ i)+ 

{ B x ,i+i )j+ i,k-\ ~ ^i+ij-i^-i)) ( A12 ) 

where B is the field on the fine grid. Note that since this is now a centered difference, the minmod 
slope limiter is not used. 
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A. 2. Implementation in Enzo 

In order to avoid complicated book keeping routines to determine which cells are being pro- 
longed into, and from which direction, we formulate only one interpolation polynomial, given by 
equations E3M The necessary finite differences for a given refinement region are taken from the 
finest data available, as in equations IA11I and IA121 The last four terms in each reconstruction 
polynomial are there exclusively to ensure consistency of Old Fine Grid Data, so for faces that have 
no Fine Data before the reconstruction, these are set to zero. Since the reconstruction polynomial 
exactly matches the old fine grid data, this also eliminates the need to copy the old fine grid data 
to the newly refined patch. 



B. Flux Correction 



At any given time in an AMR simulation, there are points in space that are described by 
more than one data structure. In a finite volume hydro calculation, with cell centered data fields, 
this occurs at the boundary between coarse and fine grids in the Flux fields, F. In an AMR MHD 
calculation, with face centered magnetic fields, this occurs at the same boundary, in the face centered 
magnetic field, and the edge centered electric field. Ensuring consistency between data is vital for 
the conservation of quantities like mass, energy, momentum, and V • B. Flux Correction is essential 
for this consistency. 



B.l. Conservation Form 

It is useful to briefly describe the basic formulation of the methods used in Enzo and EnzoMHD 
before moving on to the flux correction mechanism. 

Any conservative system, such as ideal MHD, can be written in a differential form as 

dV , . 

^ + V-F = (Bl) 

where V and F are suitably defined, in our case by [2H and [25] Here we ignore any source terms. 

In finite volume methods, we store average quantities of V and F, and re- write the conservation 
law in Conservation Form, using the Fundamental Theorem and Stokes Theorem. Starting with 
eqn IBll and integrating, we get: 



t+At f q V 

—dVdt 
v dt 



t+At 



F ■ dAdt 



(B2) 
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where the volume V is taken from the point (x, y, z) to (x + Ax, y + Ay, z + Az). Now let 

V n = ~^y J v V(x,y,z,t n )dV (B3) 
If „ 1 



2> J ' A AyAx 



F x,l+lJ,K = ~~a77~a~Z I F{x = I + -,y,z) ■ xdydz (B4) 



Aj/,A2 



where x is the unit vector in the x direction. Similar definitions apply F y and F Zl and 



F x = ir+ I F x dt (B5) 
^ -/At 



The averaging here was taken explicitly in two steps to emphasize that Ax, Ay and Az are possibly 
functions of t, as the are in cosmological hydrodynamics. Putting this all together, we get the 
equations in their final analytical form before discretization (also the last form we'll be using here) 

Vft* = Vi n ,J,K - M(^(F xJ+hJtK - F xJ _ hJtK )+ 

^ F y,i,J + hK-F yJJ . hK )+ (B6) 

Note that equation IB6I is an exact equation, since only averages and the fundamental theorem of 
calculus have been used up to this point. The trick in finite volume methods such as our MHD is 
finding appropriate approximations to F that are both accurate and stable. 



B.2. Conservation Form and AMR: Enter Flux Correction. 

As mentioned at the beginning of the section, an AMR simulation has multiple data structures 
representing a single point in space. In entirely cell centered codes such as PPM, the only such 
instance is at the surface of a fine grid boundary, where both the fine grid and coarse grid represent 
the flux at that point. Moreover, after the fine grid field is projected into the coarse, there's a 
mismatch on the coarse grid itself as to the value of the flux at the surface. The value of that 
discrepancy can be easily found. After the projection, a coarse grid at a point (I, J) has the value 
(restricting to 2d, for clarity) 

Kf= E C, +1 ( B7 ) 
i=i±\ 

3=J±\ 
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where lower case quantities denote the value of the fine grid data. Expanding the time update for 
q n+l in space and time, we find that 

_^ ^+1 _^ At" 1 - r ^~\ a At" 1 

i=I±\ m=n x,j=J±± m=n x,j=J±\ 

— (y and z terms) 

By construction of the interpolation polynomial (and projection at the last timesteps) the first term 
is just equal to V} n j, which means that, by equation IB6I Vr t effectively sees, at the point I + gj 

At ' n+1 At" 1 

-£y Fx =J2 Yl -£v^fT+l,j '=< f* > ( B9 ) 

m = n x,j=.J±\ 

However, for the cell (I — 1, J), which has no corresponding fine grid flux, F T ,i come from the 
discretization method on the coarse grid. There is absolutely no reason for the two to match, so we 
have a discrepancy in the descriptions of the data. This can be solved by simply replacing the less 
refined data that V/+i,j used with the more refined average, given by equation IB91 

At - At" 1 
Vi+i,j, /c = Vi+i,j + -gyF x j+y - Av^fZ+hi (B10) 

m j 

Now every place F x jj show up in our method, the exact same approximation is used. 



B.3. Flux Correction and MHD 

A similar formalism to that described in IB.ll is used for to advance the magnetic fields in 
EnzoMHD, but instead of using volume averages, we use area averages. The magnetic evolution is 
given by the induction equation: 

dB 

i = -" E < B11 > 

When discretized, equation IB 1 1 1 yields the equation 

At 

Ay(%, 7+ i JtK+ i -Ki+h J,K-i)) 



B n+ T 1 j = B n i — - — — (Az(E T ,i K — E j.i j i „)+ (B12) 



where 



Bl I+hjK = ^Jj(x = I+ 1 -,y,z,n.xdydz (B13) 

^ rt+At j 
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which is also exact, and the main problem is finding a suitable approximation for E. 

Again, after the area-weighted projection of the fine grid field b x into the coarse grid B x , there's 
a discrepancy between the electric field at a refined point on the sur f ace of a refined grid, as it's 



seen by both grids that have subgrids and grids that don't. In lBalsaral (|200ll ). he suggests a similar 
flux correction mechanism to that of the standard hydro, described in IB.2I However, due to an 
issue with the initial implementation of flux correction in Enzo (which has since been fixed) and 
ease of computational logic, we chose a different route. In EnzoMHD, instead of projecting fine grid 
magnetic fields into coarse magnetic fields and then correcting zones in the coarse grid, we project 
the electric field and then take the curl of the entire coarse grid. Thus, all coarse grid magnetic 
fields see the most accurate data at the same time, and no a-posteriori correction needs to be done. 
Where there are no subgrids, the coarse grid sees an electric field that comes from the CT module 
in section l2~7l and where there are subgrids it sees 

Ai 2 2 4 Z '* — 2 '3~ 2 '^+4 

While a complete flux correction treatment would potentially save on memory and flops, in practice 
the extra memory is negligible compared to the total memory and time used by the rest of Enzo, and 
the extra floating point operations done here are offset by increase cache utilization of the data, as 
the entire grid is done in a single stride one sweep instead of an essentially random access pattern. 

As described in section 12.51 some of the subgrids get their boundary conditions updated from 
the parent zones. Because of this, the curl of the magnetic field is actually taken twice. The first 
time is done immediately after the hyperbolic update, in order to ensure that the parent zones are 
up to date for the interpolation of the ghost zones of the subgrids that need it. The second time 
is after the subgrids project their electric field to the parent, to ensure maximal accuracy of the 
parent grids. This additional call takes negligible time, as the curl has relatively few operations. 
See appendix O for the details of this order of operations. 



C. Schematic for the Cosmological MHD Code 

In this section, we present a schematic of the MHD code, for clarity and easy reference. 

Step 0.- We start with conserved quantities density, total energy, and momentum (p% M , E™ otal , Pj^/), 
and primitive quantities velocity and gas pressure {v"% M , Pg as ) for the baryonic matter; face and 
cell centered magnetic fields and Lagrangian dark matter mass, position, and velocity 

(Pdm-i x ™' v Z)m)- These are all at time t n . Where needed, primitive quantities will be described 
by U = (p D M,P gas ,v D M,B), and conserved quantities by V = (pDM, EtotahPDM,^)- Conversion 
between the two is done as needed. 
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Step 1. Solve Poisson's equation for the acceleration field at t n+ 2 

r^plM+PlM (CI) 
A+n \4.n 

v 2At n ~ i; 2At n ~ 1 v 7 

n+l/2 1 fj.n+1/2 A n + 1 / 2 \ fno\ 

9i = 2a^5 Xl { ^ } (C3) 

Step 2.- Update particle positions and velocities. (Strictly speaking, this happens after the 
Expansion step, but the narrative works better if it's here.) 

\-tn -n+l/2 A+n 

n+l/2 _ n _^_a 'n At ,n+l/2 /p^ 

dm - DM 2 a^ 1 / 2 ^ ^ 

*SE = x«, M + Ai>5^V +1/2 ) (C5) 

n+ l _ n+l/2 At" d n+1 / 2 n+ l/2 _ At" n +l/2 , pR x 

V i,DM ~ V i,DM 2 a n+l/2 U i,-DM 2 Si ^ ' 

Step 3.- Apply half of the gravitational and expansion update to the fields that require it, to 
obtain the temporary state U = (p, P£, ta i, v BM , B™) 

V B M ~ V BM — — V BM — — g (L,{) 

pn = p n_^^ p n (C8) 

B™ = B?-- A B c (C9) 

U = (p,P t n otaUV BM ,B2) (CIO) 

Step 4. Compute interface states at i ± |, n + | using linear spatial reconstruction and second 
order time integration: 

U^+lv U i+}, R *= Ui-i, U U Ui+i, Ui+2 (Cll) 

Step 5. Compute approximation of the flux in equation [25] at the interface i + \. This is done 
by solving the Riemann problem using one of the solvers mentioned in section 12.61 

F. \ = Riemann(U. \ 2 T , U. \\) (C12) 

^+2 i- ^^2' ^+2;-^ 



Step 6. Update the conserved quantities with the new fluxes: 

At 



(r)M flD = - ^[F i+k - (C13) 
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Step 7. Compute Electric field from Fluxes 



i 



<=*i+i (C14) 

Step 9. Update magnetic fields from electric fields for the first time. 

B n f +l = B n f - — V x E n ^ ! (C15) 

Step 8. -Gravitational step for the baryonic matter, with time centered density 

( n n -I- n n+1 ) 

i n+1 \ _ n _ A.n \^_^J\ihd]_ n+1/2 / nifi x 

\Vi,BM>MHD,Grav ~ \*?i,BM ) MHD LXl 2 ^ ^ 1U 7 

Step 9.-Expansion step for the baryonic matter, 

(v n. + n = 1 - (AtV2)(a"+VVa^) +1 

V V BMJMHD,Grav,exp j + ^Vq/I+I^ \ V BM >MHDGrav \^ 1 1 ) 

x i_(A f ")(d"-n/y ^) 

P l + (At»)(d"+V2/ a n+l/2)^ ) mhd 

Step 10. Recurse to finer grids. Integrate fine grids from t n to 

^F-LeGrids < < = ^FineGrids (C19) 

Step 11. -Flux correction step for conserved baryon field quantities 

T/-n+l - /pn+l/2\ /pn+l/2\ T/ n + 1 fPOD"! 

%lifD,ft«-,. : xp,/c ^ ^ /' JFineGnds, V MHDG rav,exp {^ZV) 

Step 12.-Project conserved baryon field quantities and electric field from fine grids to coarse 
grids. This is done after the flux correction step to avoid any bookkeeping errors. The average is 
taken over At™ and the surface of each FineGrid. 

VparentGrid = < ^FineGrid >t,surface (C21) 
EparentGrid = < ^FineGrid >t,surface (C22) 



Step 13. Update magnetic fields from electric fields for the final time 

At 1 
a 



B] +1 = B] - x E P J entGrid (C23) 



Step 14. Apply expansion to the Face Centered Fields 



?n+1 _ 1 - (AP/4)(«" +1 / 2 /« n+1/2 ) , „n + i, 



i + (M n /A){a n + l l 2 /a n+l / 2 y f ' [ ' 
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Step 15. Compute cell centered magnetic field from face centered (with the expansion subscript 
from step 9 dropped for clarity) 

B ciXj,k = °- 5 * ( B f,x,t+ij,k + B f, x ,i-i,j,k) 

Step 16. We have now finished an update of this level. Rebuild the hierarchy from this level 
down. 

VnIw FineGrids ^=V n+l (C26) 
B f,New FineGrids < ^ = ^\f (C27) 
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